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Normal-conducting mesoscopic systems in contact with a superconductor are classified by the 
symmetry operations of time reversal and rotation of the electron's spin. Four symmetry classes are 
identified, which correspond to Cartan's symmetric spaces of type C, CI, D, and DHL A detailed 
study is made of the systems where the phase shift due to Andreev reflection averages to zero along 
a typical semiclassical single-electron trajectory. Such systems are particularly interesting because 
they do not have a genuine excitation gap but support quasiparticle states close to the chemical 
potential. Disorder or dynamically generated chaos mixes the states and produces novel forms of 
universal level statistics. For two of the four universality classes, the n-level correlation functions 
are calculated by the mapping on a free ID Fermi gas with a boundary. The remaining two classes 
are related to the Laguerre orthogonal and symplectic random-matrix ensembles. For a quantum 
dot with an NS-geometry, the weak localization correction to the conductance is calculated as a 
function of sticking probability and two perturbations breaking time-reversal symmetry and spin- 
rotation invariance. The universal conductance fluctuations are computed from a maximum-entropy 
S-matrix ensemble. They are larger by a factor of two than what is naively expected from the analogy 
with normal-conducting systems. This enhancement is explained by the doubling of the number of 
slow modes: owing to the coupling of particles and holes by the proximity to the superconductor, 
every cooperon and diffuson mode in the advanced-retarded channel entails a corresponding mode 
in the advanced- advanced (or retarded-retarded) channel. 
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I. INTRODUCTION 

Following the early work of Wigner jl), Dyson in his 
classic 1962 paper [0 classified complex many-body sys- 
tems such as atomic nuclei according to their fundamen- 
tal symmetries. Arguing on mathematical grounds, he 
proposed the existence of three symmetry classes, which 
are distinguished by their behavior under reversal of the 
time direction and by their spin. The statistical proper- 
ties of these classes arc described by three random-matrix 
models, called the Gaussian Orthogonal, Unitary and 
Symplectic Ensembles (GOE, GUE, GSE). Dyson's clas- 
sification scheme has since proved very far reaching. Al- 
though atomic nuclei display only GOE statistics, physi- 
cal realizations of the other two classes were later found in 
chaotic and disordered single-electron systems subject to 
a magnetic field (GUE) or to spin-orbit scattering (GSE). 

By standard arguments, Wigner-Dyson statistics ap- 
plies to the ergodic limit, i.e. to times long enough for 
the degrees of freedom to equilibrate and fill the avail- 
able phase space uniformly. More specifically, in the con- 
text of disordered mesoscopic systems the ergodic limit 
is reached for times larger than the diffusion time L 2 /D, 
where D is the diffusion constant and L the linear ex- 
tension of the system. By the uncertainty relation, the 
ergodic limit corresponds to the energy range below the 
Thouless energy fiD/L 2 . 

One may ask whether the level statistics of disordered 
or chaotic single-particle systems in the ergodic limit 
must always be Wigner-Dyson or whether different statis- 



tics is possible. The answer is that Wigner-Dyson statis- 
tics is generic and universal as long as the statistics is re- 
quired to be stationary under shifts of the energy. (This 
can be understood from the mapping on a nonlinear a 
model ||.) However, if the stationarity condition is re- 
laxed and additional symmetries are imposed, new uni- 
versality classes may arise. This happens, for instance, 
when a massless Dirac particle is placed in a random 
gauge field. Because the Dirac operator D anticommutes 
with 75 in the chiral (or massless) limit, its matrix is 
block off-diagonal in the eigenbasis of 75. As a result, the 
eigenvalues of D are either zero or come in pairs (A, — A). 
The average spectral density of D close to zero is nonsta- 
tionary but universal and is of relevance for the physics 
of QCD at low energies. It is determined by one of three 
so-called chiral Gaussian ensembles, where different en- 
sembles correspond to different choices of the gauge group 
and the number of flavors jij . These ensembles have ap- 
peared in the context of disordered single-electron sys- 
tems, too ||. 

In the present paper we introduce and analyse four ad- 
ditional Gaussian random-matrix ensembles, which share 
many striking similarities with the chiral ones but are 
demonstrably distinct. The universality classes they de- 
scribe are realized in mesoscopic NS-systems, i.e. in 
microstructures composed of a metallic part in contact 
with one or several superconducting regions. Just as in 
the classic Wigner-Dyson case, the universality classes 
are distinguished by their behavior under time rever- 
sal and rotation of the (electron's) spin. The four new 
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classes together with the six known ones add up to a 
grand total of ten. We have reasons to believe that this 
exhausts the number of possible universality classes in 
disordered single-particle systems and none else will be 
found. (More precisely speaking, by universality classes 
we here mean infrared RG fixed points describing an er- 
godic limit.) Some of our ideas were anticipated in ||[7]]. 

The prototype of the kind of system we are going to 
study is depicted in Fig. |l|. A metallic (i.e. normal- 
conducting) quantum dot is put in contact, via potential 
barriers, with two superconducting regions. Several leads 
are attached for the purpose of making current and volt- 
age measurements. The metallic quantum dot may or 
may not be disordered. In the latter case we assume its 
geometric shape to be such that the classical motion of 
a single electron inside it is chaotic. The quantum dot 
may be pierced by a magnetic flux of the order of one or 
several flux quanta, and there may exist some impurity 
atoms causing spin-orbit scattering. The temperature is 
so low that the electron's phase coherence length exceeds 
the size of the quantum dot by far. 




FIG. 1. Metallic quantum dot (N) in contact with two su- 
perconducting regions (S). The dot is separated from the leads 
(L) by a tunnel barrier (T). 

The characteristic feature that distinguishes the above 
kind of quantum dot from more conventional mesoscopic 
systems, is the possibility for two electrons to tunnel 
through the potential barrier at the NS-interface, thereby 
adding a Cooper pair to (or removing it from) the su- 
perconducting condensate. An equivalent statement in 
single-particle language is that an electron incident on 
the NS-interface may be retroreflected as a hole (and 
vice versa). This process of particle- hole conversion, 
which conserves energy, momentum and spin but vio- 
lates charge, is called Andreev reflection. In the semi- 
classical limit, Andreev reflections give rise to numerous 
almost-periodic orbits whose action does not grow but 
remains of order h as the length of the orbit increases 
H . The existence of these orbits modifies the mean den- 
sity of states (Weyl term) of the quantum dot without 
leads: in general, an excitation gap opens and we arrive 
at the "boring" situation where the vicinity of the chemi- 
cal potential is devoid of single-particle states. However, 



by tuning the phase difference of the order parameters 
of the two superconducting regions to the special value 
7r, we can make the gap close. More generally, we expect 
quasiparticle excitations to exist right at the chemical po- 
tential whenever the phase shift incurred during Andreev 
reflection vanishes on average over the NS-interfacial re- 
gion. Disorder or dynamically generated chaos mixes the 
states and creates a universal spectral region close to the 
chemical potential. Its width is determined by the energy 
uncertainty which is caused by the coupling of particles 
and holes by Andreev reflection. It is this very region and 
its consequences for the transport properties that we are 
going to study in the present paper. 

The organization of the paper is as follows. Mesoscopic 
independent-quasiparticle systems are classified accord- 
ing to their behavior under time reversal and spin rota- 
tions in Sec. pi H aving specified the required dynamical 
input in Sec.^IIl| we formulate the appropriate random- 
matrix ensembles in Sec. [fv|. In Sec. we discuss the 
spectral statistics of an isolated system, using first the 
Dyson-Mehta orthogonal polynomial method and then 
diagrammatic perturbation theory. The latter method 
easily extends to the calculation o f th e transport prop- 
erties of an open system. In Sec. VII we work out the 
weak localization correction to the average conductance 
and in Sec. VIII the universal conductance fluctuations. 



Our conclusions are presented in Sec. 



II. SYMMETRY CLASSIFICATION 

The treatment of this paper is based on the BCS 
Hamiltonian in the Hartree-Fock-Boboliubov mean-field 
approximation: 

H= fd d xl Y, + A^Vl + A*ViV>t I , 

\o-,t=T4 / 
h = (p - eA) 2 /2m + V + U so ■ <x x (p - e A) - /i. 

Here V(x) is a scalar potential which may have a random 
component, and A(x) is the pairing field. The presence 
of a magnetic vector potential A(x) breaks time reversal 
symmetry while the spin-orbit field Uso(x) breaks in- 
variance under rotations of the electron's spin. \i is the 
chemical potential. 

The second-quantized Hamiltonian H can be rewritten 
in an equivalent first-quantized form by the Boboliubov- 
deGennes (BdG) formalism. For our purposes it is con- 
venient to introduce some generic orthonormal basis of 
single-electron states \a), where a is a multi-index that 
combines the orbital and spin quantum numbers of the 
electron. If N is the number of orbital states used, a runs 
from 1 to 27V. Let c' a and eg be the usual creation and 
annihilation operators obeying the canonical anticommu- 
tation relations cj,c^ + c^c\ = 5 a p- The Hamiltonian H 
can be written 



2 



a/3 



Hermiticity requires h a p — h% a , and the matrix ele- 
ments A a p must be antisymmetric by Fermi statistics: 
A a f} = —Af} a . Now we write H in the form "row multi- 
plies matrix multiplies column" : 

H = ~ ( c t c) ( \* ) ( c C t ) + const. (1) 

In this way every Hamiltonian H is uniquely assigned to 
a AN x 4iV-matrix TL, 



n 



h 

-A* 



A 

-h T 



(2) 



The eigenvalue problem for TL is known as the 
Bogoliubov-deGennes equations. We refer to TL as the 
"BdG-Hamiltonian" for short. 

The first-quantized Hamiltonian TL acts in an enlarged 
space, namely the tensor product of the physical space 
C 2N (orbitals and spin) with an extra degree of free- 
dom C 2 , which we call the "particle-hole space". Note 
however that the "particles" and "holes" of the BdG- 
formalism are not the particle and hole states of a de- 
generate Fermi gas. Indeed, the matrix h already acts 
on all of the single-electron states, which have energies 
either above or below the chemical potential. The BdG- 
hole states acted upon by — h T are identical (and in 
this sense redundant, or unphysical) copies of the BdG- 
particle states acted upon by h. They are introduced for 
the convenience of treating the pairing field within the 
formalism of first quantization. 

The aim of the current section is to classify systems ac- 
cording to their symmetries. Using the BdG-formalism 
we will show that the presence or absence of time-reversal 
and/or spin-rotation invariance leads to four distinct 
symmetry classes. The situation thus is different from 
the well-known Wigner-Dyson scenario where only three 
distinct classes exist. 

The discussion of Sees. 



IIA-IID uses some basic facts 



from the theory of Lie algebras and symmetric spaces 
and is somewhat technical. The casual reader may wish 
to skip these subsections and proceed directly to Table 
1 given at the end of Sec. [ID . A brief summary of the 
symmetries of TL for each class is provided also at the 
beginning of Sec. [V. 



skew-symmetry of A. Because the set of hermitian ma- 
trices does not close under commutation whereas the an- 
tihermitian ones do, we prefer to work with X := iTL 
rather than TL in the current section. In terms of A", the 
conditions h = h\ A = — A T can be presented summar- 
ily in the form 



X T = X 







-Yj x X Yjfr 

1-2AT 




(3) 



"-2JV- 



If two matrices X, Y satisfy these equations, then so does 
their commutator [X, Y]. Hence, we may view X as an 
element of some Lie algebra. To identify this Lie alge- 
bra we conjugate by X i— > X = UqXUq 1 where Uq — 

"73 ( \ \ 1 ® 1-2/v- Equations (||) then take the form 

= X = —X T or, equivalently, X = X* = -X T . 
This shows that (Q) fixes a so(47V)-algebra, i.e. a Lie 
algebra isomorphic to the real antisymmetric AN x AN- 
matrices. Since so^A^) = Z?2/v in Cartan's notation, we 
denote the present symmetry class by the symbol "D" . 

Being a Lie algebra element, X can be brought into 
diagonal form by I h> 1! = gXg^ 1 where g is an ele- 
ment of the corresponding Lie group which is isomorphic 
to SO(4A0 and is defined by 



-it ^ -i T ^ 

; =9 = Z x g S 2 



(4) 



The conditions ([}]) imply f2 = diag (iuj, —iuj) — a z ® iuo 
where lj = diag(ui, LU2, ...,L02n) with real Wj. The condi- 
tions for the canonical anticommutation relations to be 
invariant under a transformation 



-t 



7 



7 



can be shown to coincide with (Q). Thus, inserting 
X = g~ x £lg into (lj) we obtain 



H 



w a I7a7a -7 A 7l 



with 7^7^ + 7^7^ = S a p- The frequencies u>\ may be 
positive or negative. The BCS ground state is defined by 
demanding 7a|BCS) = for uj x > and 7l|BCS) = for 
lo\ < 0. The normal-ordered Hamiltonian 



H : = 



J2 i^it1ta+ i^itat! 



w A >0 



"A<0 



A. Symmetry Class D 



is always positive. 



We start by considering systems with the least degree 
of symmetry, i.e. systems with neither time-reversal sym- 
metry nor spin-rotation invariance. In this case the ma- 
trices h and A in general have no symmetry properties 
beyond those stated above, namely hermiticity of h and 



B. Symmetry Class C 

We now consider systems without time-reversal sym- 
metry but with spin-rotation invariance. We again 
use the unique representation of a second-quantized 
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BCS mean-field Hamiltonian H by (i times) a BdG- 
Hamiltonian X = iTL. 

We write the particle-hole decomposition of A as A = 

^ ^ ^ ^ or, in tensor-product notation, 

X = E pp <g> A + E ph <g> B + E hp <g> C + E hh <g> D. 

The condition X = -£ X A T £ X means B = -B T , C = 
— C T and D = — A T . Antihermiticity requires A = —A'' 
andC^-Bt. 

The generators of spin rotations, Jk (k = x,y,z), are 
represented on particle-hole space by Jk — (E pp <£> Cfc — 
.Ehh <8> cj) ® ljv- Spin-rotation invariance of the Hamil- 
tonian requires that X and Jk commute. This condi- 
tion is easily seen to constrain A, B, C to be of the form 
A = 1 2 <8> a, B = i<T y ® b and C = —%o y ® c or, in matrix 
presentation, 



(a 








6 ^ 





a 


-b 








—c 


T 

-a 1 





Vc 








-a T 7 



We see that X decomposes into two commuting sub- 
blocks. One corresponds to spin-up particles and spin- 
down holes, and the other to spin-down particles and 
spin-up holes. Because the subblocks are related by an 
algebra homomorphism (b —6, c — c) it is sufficient 
to focus on one of them, say 

X ^=( a c -a T ) 

and account for spin degeneracy by inserting factors of 
two whenever needed. We drop the subscript and write 
AforA T . 

Since B = — £> T , the equation B = ia y ® b implies 
b = +b T . Similarly, we deduce c = +c T . Antihermiticity 
requires a — — and c = — 6L All these conditions are 
summarized by 

- A* = X = -E y X T E y , Y, y = a y ®l N . (5) 

This is the defining equation of the symplectic Lie al- 
gebra sp(2A). Thus X = iH is an element of sp(2A). 
Since sp(2A) = Cn in Cartan's notation, we denote the 
present symmetry class by "C" . 

The second-quantized Hamiltonian associated with 

* = Me-« T )' S 

2 ^ V Cmn anm J \ C nl/ 

As before, we can diagonalize X by X = g^ x ilg where 
fi = <J Z ® and oj — diag(wi, ...,u>m), and g now sat- 
isfies 5 ~ lf = ,g = £ y . 9 - lT £ y , i.e. 5 e Sp(2A). The 
transformation 




diagonalizes the Hamiltonian: 
if = _ff | + 

= ^ Z WA (Ja T 7a T + 71i7ai - 7a|7a| - 7 Ai 7A|) • 

A 

Because g £ Sp(2A) C U(2A) is a unitary matrix, the 
canonical anticommutation relations are preserved by the 
transformation from (c,c^) to (7,7^). Every quasiparti- 
cle level has a trivial degeneracy due to spin. 

C. Symmetry Class Dili 

We next consider systems with time-reversal symme- 
try but without spin-rotation invariance. Recall that 
the conditions for symmetry class D, — X^ = X = 
-£ X X T £ X with E x = a x <g> 1 2 ® ljv, fix a so(4A)- 
algebra. The time-reversal operator T acts on the BdG- 
Hamiltonian by 

where r = l^®^® ljv- Using X = iH and X* = — X T 
we get the action of T on A, T : X tX t t~ 1 . Thus, 
for a time-reversal invariant system, X is subject to the 
additional constraint A = +tX t t~ 1 . We denote the set 
of solutions in so(4A) of this condition by V . While V 
does not close under commutation and therefore does not 
form a subalgebra of so(4A), the solution set, JC, of the 
complementary condition Y = —tY^t^ 1 does. There- 
fore, we may describe V as the complement of a subalge- 
bra JC in so(4A). In formulas, so(4A) = JC + V. We are 
now going to identify JC. 

The equations for JC can be rewritten 

-Y^=Y= -^ X Y T ^ X = -(S x r)F(S x r)- 1 . 

Let Uq be the unitary matrix given in particle-hole and 
spin decomposition by 

Conjugation by [/ , FhF= Uq 1 YUq, takes the equa- 
tions for IC into 

== ^Y == 5j x ^^ = ^j^YYiz. 

The solutions of the latter are of the form Y = 
diag (Z, —Z T ^j with Z an antihermitian 2 A x 2A-matrix. 
We now recognize JC as being isomorphic to the Lie alge- 
bra of antihermitian 2 A x 2A-matrices, or JC ~ u(2A). 
Thus, the space V of BdG-Hamiltonians A = iH is ob- 
tained from so(4A) by removing a u(2A)-subalgebra. 
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Because V is the difference of two Lie algebras so(4Af) 
and u(2iV), it can be interpreted as the tangent space 
of the quotient SO(47V)/U(27V) of the corresponding Lie 
groups, which is a symmetric space of type Dili in Car- 
tan's notation; hence the name Dili for the present sym- 
metry class. 

From the general theory of symmetric spaces we 
know that an element X G V can be diagonalized by a 
transformation X i— > k~ 1 Xk with k £ exp/C = U(2iV). 
Time-reversal symmetry causes every eigenvalue to be 
doubly degenerate by Kramers' theorem. 



D. Symmetry Class CI 

Finally, we turn to systems with both time-reversal 
symmetry and spin-rotation invariance. Recall that spin- 
rotation invariance causes the BdG-Hamiltonian X 
1 1 1 ■ to obey the relations — X' = X = — T, y X T 'Sy, 
which define the symplectic Lie algebra sp(2A^). Because 
of the restriction to a single spin sector, the action of 
the time reversal operator simplifies to T : I h X T . 
Thus, time-reversal symmetry constrains X to be sym- 
metric. Let K, now denote the subalgebra of antisymmet- 
ric matrices in sp(2iV). Then X, being symmetric, lies 
in the complementary set V defined by sp(2A r ) = JC + V . 
We claim that K, is isomorphic to the unitary Lie alge- 
bra u(iV). To prove this, we observe that the solutions 
Y e K. of -yt = y = -T, y Y T E y = y T have the form 
I2 ® Re^4 + i<jy ® IvaA where A is an arbitrary antihcr- 
mitian N x iV-matrix, i.e. A S u(iV). The embedding 
u(N) i-> sp(2iV) by A i-> 1 2 ® ReA + ia v ® ImA is eas- 
ily seen to be a homomorphism of Lie algebras. There- 
fore K. ~ u(N) as claimed. The linear complement V of 
u(N) in sp(2A r ) can be regarded as the tangent space of 
Sp(27V)/U(A r ), which is a compact symmetric space of 
type CI according to Cartan's list. 

For the benefit of the casual reader the various sym- 
metry classes and the names by which they are referred 
to in the present paper, are summarized in Table 1. 



Class 


time-rev. 


spin-rot. 


sym. space 


D 


no 


no 


SO(4A0 


C 


no 


yes 


Sp(2iV) 


Din 


yes 


no 


SO(4A0/U(2iV) 


CI 


yes 


yes 


Sp(2iV)/U(A0 



Table 1 



E. Is the symmetry (g|) wiped out by Coulomb 
effects? 

The symmetry (^) is central to our approach. Just how 
robust is it? 

The relations (||) follow from the well-known mathe- 
matical fact flOl that the set of all bilinear combinations 



of the fermion creation and annilihation operators is iso- 
morphic to an orthogonal Lie algebra. Put differently, 
the symmetry (^|) requires no more than the fermionic 
nature of the electron and the use of the Hartree-Fock- 
Bogoliubov mean-field approximation, allowing us to ex- 
press the Hamiltonian in terms of bilinears of the creation 
and annihilation operators. Alternatively, we could say 
that (||) is an exact symmetry whenever the system can 
be described in terms of independent Boboliubov quasi- 
particles. 

What happens when we add a Coulomb charging en- 
ergy to the Hamiltonian? The relative minus sign be- 
tween the particle-particle and hole-hole blocks of H, 
Eq. (Q), tells us that, if the creation of an electron in 
a given state costs a certain amount of energy, then the 
creation of a hole (removal of an electron) in this state 
should release exactly the same amount. The Coulomb 
interaction, however, does not conform to that princi- 
ple. When a charge is added to a charge-neutral sys- 
tem, say, it makes no difference whether this charge is 
a particle or a hole, the electrostatic energy cost is pos- 
itive in both cases. Therefore, the Coulomb charging 
energy (as well as other perturbations that do not fit 
into the independent-quasiparticle framework) violates 
the symmetry (|^). More precisely speaking, we expect 
the independent-quasiparticle approximation to be ade- 
quate for describing the short-time physics, but at suf- 
ficiently long times Coulomb effects must become vis- 
ible and, in particular, they will cut off the particle- 
hole modes we are going to study in the present pa- 
per. Whether the cutoff time can be long enough for the 
consequences of the symmetry (||) to be observable, is a 
tough quantitative question for theory, which cannot be 
answered without an understanding of screening in open 
and finite metallic systems. Fortunately, the question 
has already been answered in the affirmative by experi- 
ment. Over the last years a number of novel mesoscopic 
NS-phenomena has been observed, the most prominent 
of which is the dramatic enhancement [Q of the differ- 
ential conductance at zero bias in NS-geometries with a 
high potential barrier separating the normal-conducting 
and superconducting regions. This phenomenon has been 
explained p"3[ ] by a mechanism called "coherent Andreev 
reflection" or "reflectionless tunneling" juj , which is the 
result of constructive interference between semiclassical 
paths with one Andreev reflection and a variable number 
of normal reflections. In order for such an interference 
to take place, the dynamical phases of a particle and 
a hole traversing the same path in opposite directions 
must cancel each other. It is precisely the symmetry (|j|) 
in combination with the extra symmetries defining class 
CI, that guarantee the necessary phase relation between 
particles and holes to hold. We conclude that there exists 
convincing experimental evidence that the symmetry (|^) 
is not wiped out by the Coulomb interaction but leads to 
observable consequences. In the remainder of this paper 
we are going to ignore Coulomb effects. 
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III. DYNAMICAL INPUT 

The classification of Sec. [H] refers only to symmetry 
and thus is very general. To go further, we make two 
dynamical assumptions. 

When an electron is retroreflected from the NS- 
interface as a hole, its wavefunction acquires a phase shift 
which is determined by the phase of the superconducting 
order parameter. Our first assumption is that this phase 
shift, here called the "Andreev phase shift" for short, 
vanishes on average over the NS-interfacial region. To 
appreciate what such an assumption implies, let us look 
at a few examples. Consider first an SNS-system con- 
sisting of an infinite slab of normal metal sandwiched 
between two superconducting slabs Si and £2- The pair- 
ing interaction causes the existence of an excitation gap 
in each of the superconducting regions. We now ask how 
the presence of the normal-conducting slab affects the 
excitation spectrum of the combined SNS-system at the 
chemical pot ential. The answer to this question was first 
given in |[15| , fL6[ and it is essentially as follows. In the 
clean limit, the BdG-Hamiltonian is separable and we 
can get a qualitative understanding of the quantal en- 
ergy spectrum by the method of semiclassical quantiza- 
tion. For simplicity we assume that all reflections at the 
NS-interface are Andreev. Every periodic classical mo- 
tion then is some multiple of a primitive periodic orbit 
where an electron moves back and forth between the su- 
perconducting regions and is Andreev reflected at each 
interface. If k p (fch) are the wave numbers of the par- 
ticle (hole) normal to the slabs and L is the thickness 
of the normal-conducting region, the Bohr-Sommerfeld 
quantization rule applied to this periodic motion reads 

± (fc p — kh)L + 7r + ip 1 — (f2 = 2ttu (n 6 Z). (6) 

Here it + fi — f2 is the phase accumulated by the two 
Andreev reflections if f\ and f2 are the phases of the su- 
perconducting order parameter in the regions S± and S2 ■ 
For an electron with energy equal to the Fermi energy, k p 
equals fch, so the first term on the left-hand side vanishes. 
From the resulting equation f\~f2 = (2n — 1)tt we see 
that the quantization condition can be fulfilled only when 
f 1 and f 2 differ by an odd multiple of tt. In other words, 
for cos((/9i — V2) 7^ — 1) which includes the homogeneous 
case if 1 —f2, we expect a gap in the excitation spectrum 
not only in the superconductor but also in the combined 
SNS-system. On the other hand, for cos((^i — ^2) — > — 1 
the gap closes and quasiparticle excitations exist all the 
way down to zero energy. The latter situation is special 
in that the Andreev phase shift vanishes on average over 
the two NS-interfaces in that case. 

The above argument applies to the extreme limit of a 
clean system which clearly is unrealistic. What can we 
say about the effects of disorder? A generic random po- 
tential destroys separability and makes Bohr-Sommerfeld 
quantization inapplicable. The general case therefore 
needs to be studied with the help of a computer, or by 



using the random-matrix theory that will be developed 
in the remainder of our paper. What is easy to treat ana- 
lytically is the case of a slowly varying random potential 
depending only on the coordinate, z, of the direction per- 
pendicular to the slabs. In this case the quantization rule 
(H) remains valid if we replace the expression (fc p — fch)Z/ 

by the integral f$ (kp(z) ~ k^(z)jdz where k p (z) — 

pmGu + E-FO*)] 1 / 2 , fch(z) = [2m(^-£-V A (z)] 1 /2 ; an d 
E = \i + e is the total energy of the electron. Since 
fcp(z) = fch(z) for e — 0, our conclusions are the same as 
before: there is a gap for cos(y>i — if 2) 7^ — 1, and the gap 
closes as cos(y>i — f>2) — * — I- 

Another instructive example is provided by the vortex 
solution for a clean type-II superconductor. The phase 
of the superconducting order parameter uniformly winds 
once around the unit circle as we go once around the 
vortex core. For this reason, the pairing field experi- 
enced by normal excitations bound to the vortex core 
vanishes on average over the vortex. Because the vortex 
solution breaks translational symmetry, there must ex- 
ist some RPA (or vibrational) zero modes of the vortex. 
These zero modes are the Goldstone modes associated 
with the spontaneous breaking of translational symme- 
try by the localized vortex solution. It follows that, if 
the RPA correlation energy vanishes (is small), i.e. if 
the excitation energies are given by (are approximately 
given by) sums of two quasiparticle energies, there must 
exist quasiparticle excitations with vanishing (small) en- 
ergy. In contrast, for a piece of cylindrically shaped nor- 
mal metal immersed in a superconducting environment 
("columnar defect") there is no general reason why we 
should expect quasiparticle excitations with low energy. 

These two examples, the SNS slab geometry and the 
vortex, lend support to the plausible expectation that 
a pairing field which is locally nonzero but whose mean 
phase (e 1 ^) vanishes in a suitably defined sense, is inef- 
fective at creating a genuine gap in the density of states 
near the chemical potential. This then is the motivation 
behind the above requirement that the Andreev phase 
shift should vanish on average over the NS-interfacial re- 
gion: it ensures that the gap closes and quasiparticles 
can exist right at the chemical potential. 

Our second main input is the assumption that the clas- 
sical dynamics in the normal- conducting (N-) region be 
chaotic. The presence of a sufficient amount of disor- 
der will always guarantee this assumption to be justi- 
fied. For a ballistic system, chaotic dynamics is achieved 
by choosing for the boundary of the N-region some sur- 
face that causes a typical classical trajectory to be un- 
stable. Chaoticity of the classical motion means that the 
long-time behavior of the system is unpredictable; in par- 
ticular, the phase shifts acquired by Andreev reflection 
along a typical semiclassical trajectory form a random 
sequence. This randomness will allow us to model the 
pairing field by a stochastic variable. Note that the spa- 
tial constancy of the magnitude of the pairing field in the 
bulk of the superconductor is an irrelevant feature for our 
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purposes. If we switch from coordinate representation to 
a generic basis of single-particle states, say the eigenbases 
of h and — h T , both the phase and the magnitude of A 
will fluctuate and be distributed around zero. 

Consider now an isolated finite system, so that the 
Boboliubov quasiparticle spectrum is discrete. According 
to our above arguments, we expect the existence of levels 
close to the chemical potential in the pure system under 
the conditions described. The effect of dynamically gen- 
erated chaos and/or disorder will be to cause mixing of 
these levels. For the conventional N-system, such mixing 
is known to lead to universal level statistics, depending 
only on symmetry. (More precisely, the level correlations 
are universal in the low-frequency regime corresponding 
to the long-time or ergodic limit.) For the case of disor- 
dered metallic granules, the level correlations were calcu- 
lated by Efetov || . His results show that the level statis- 
tics is Wigner-Dyson, i.e. identical to that of an ensem- 
ble of random matrices with uncorrelated Gaussian dis- 
tributed matrix elements. In the NS-systems considered 
in the present paper, new features appear at low excita- 
tion energy, owing to the coupling of particles and holes 
by the process of Andreev reflection at the NS-interface. 
As was shown in Sec. JO], the presence of the pairing field 
A leads to novel symmetries. We therefore expect new 
types of universal level statistics to emerge in such sys- 
tems. The new type of correlation will extend over an 
energy range set by the energy uncertainty due to the 
action of the pairing field (or Andreev reflection). The 
goal of our paper is to give a quantitative description of 
precisely these correlations and their effect on the trans- 
port properties. To reach this goal we may follow two 
different routes. The first and more comprehensive one is 
to generalize Efetov's analysis, i.e. to construct an effec- 
tive field theory of the nonlinear a model type and solve 
the field theory in the zero-dimensional limit correspond- 
ing to the universal regime. Such an approach yields not 
only the universal behavior but also the crossover to the 
short-time regimes. Since our interest is in the univer- 
sal limit, there exists also another option. Armed by the 
experience gained from the study of the N-system, we 
may replace the BdG-Hamiltonian (|2|) by an ensemble of 
random matrices with maximum entropy, paying atten- 
tion only to the symmetries under time reversal and spin 
rotation. While the field-theoretic method is more versa- 
tile, the random-matrix or maximum-entropy approach 
has the great advantage of being much simpler techni- 
cally For this reason we have chosen to follow the latter 
in the present paper. 

To maximize the entropy of the random-matrix en- 
semble we will take the matrix elements of the BdG- 
Hamiltonian to be normally distributed and statistically 
independent. All matrix elements will be chosen to have 
zero mean. For the off-diagonal blocks of the BdG- 
Hamiltonian, this choice corresponds to our assumption 
that the spatial average of the Andreev phase shift van- 
ishes on average. In general, we would need to distinguish 
between the strength of fluctuation of h and A. However, 



at low energy, i.e. within the energy window defined by 
the uncertainty due to Andreev reflection, this distinction 
turns out to be irrelevant and we may take the strengths 
to be equal. The resulting random-matrix ensemble de- 
pends on two parameters only. These are the strengths 
of the perturbations that break time-reversal and spin- 
rotation invariance and are responsible for the crossover 
between universality classes. 

To finish off this orientational section, we mention an- 
other realm of application of the above random-matrix 
ideas. Consider an array of superconducting grains or 
islands embedded in a metallic (non-superconducting) 
host. The grains are disordered and/or of irregular shape, 
and they are mutually coupled by Josephson tunnelling. 
The array is exposed to a subcritical magnetic field which 
penetrates the host but is ejected from the grains. By 
tuning the strength of the field we can frustrate the cou- 
pling between the grains and drive the system into a spin- 
glass type of phase where superconducting order exists 
locally but not globally. Such a system has been called 
a superconducting glass 0. Its prime characteristic is 
that the pairing field, or superconducting order parame- 
ter, continues to be nonzero on each grain but vanishes 
on average over large scales. The low-energy quasipar- 
ticle excitations of such a system are predicted to be 
described by the random-matrix model formulated be- 
low. Because of the breaking of time-reversal symmetry 
by the magnetic field, the relevant symmetry class is C. 
The presence of spin-orbit interactions causes crossover 
to D. 



IV. RANDOM-MATRIX ENSEMBLES 

To prepare the formulation of the random-matrix en- 
sembles, we summarize the discussion of Sec. || by pre- 
senting the symmetries of the BdG-Hamiltonian Ti for 
each symmetry class explicitly. 

For systems where all symmetries are broken (class D) 
Ti satisfies Ti = — Y, x Ti T Y, x with S x = a x ®\2®^-N ■ The 
block decomposition 



n 



A 
St 



B 

-A T 



expresses the particle-hole structure of Ti. The off- 
diagonal block B is antisymmetric by Fermi statistics. 
Hermiticity of the Hamiltonian requires A = . 

Class Dili consists of the systems where time reversal 
is the only good symmetry. For such systems Ti obeys the 
additional relation Ti — tH t t~ 1 with r = I2 <£> i<r y ® In- 
The decomposition of Ti according to particles and holes 
(outer block structure) and spin (inner block structure) 
has the form 



Ti = 



a it 
b 1l 



a U 

T 

A 

— Ot-f 



6 TT 


b u 




b ll 


T 


a n 






— Cbf-r 



7 



with antisymmetric a-^, a^, fejj and b^. Hermiticity 



,t 

Z -LT' 



and 



of 7i requires a-ff 

b n = - b \r 

For class C spin is conserved while time-reversal sym- 
metry is broken. In this case TL commutes with the spin- 



al 



hh 



N 



rotation generators Jk = (E pp eg) <Jk 
or, equivalently, TL obeys TL = JkTLJk- The particle-hole 
and spin decomposition of TL reads 



TL 



a 











a 







-6t 


T 

—a 


&t 









with symmetric b. Every level has a trivial twofold degen- 
eracy due to spin. Without loss of information we may 
focus on the spin-up sector with reduced Hamiltonian 



= ot. 



Hermiticity requires a 

In class CI both spin rotations and time reversal are 
good symmetries. The BdG-Hamiltonian satisfies TL — 
tTL t t~ 1 — JkTLJk, and is constrained by these symme- 
tries to be of the form 



TL 



fa 











a 


-b 





-b 


T 

—a 


\b 









with symmetric a and b. Hermiticity then implies that a 
and b are real matrices. 

Now recall the dynamical conditions formulated in 



Sec. III. By assumption, the classical dynamics in the N- 
system is chaotic and the Andreev phase shift vanishes on 
average over the NS-interfacial region. We therefore may 
replace the BdG-Hamiltonian by a random matrix (of the 
appropriate symmetry) with matrix elements that have 
zero mean. The principle of least information, or maxi- 
mum entropy, then leads us to postulate a random-matrix 
ensemble with a Gaussian probability distribution 



cxp (-TrTL 2 /2v 2 ) dTL 



(7) 



for each symmetry class. Here dTL denotes a Euclidean 
measure on the linear space of BdG-Hamiltonians with 
metric Ti(dTL) 2 . 

More generally, we can formulate a two-parameter fam- 
ily of Gaussian random-matrix ensembles which interpo- 
lates between all four symmetry classes. Because a Gaus- 
sian distribution (with zero mean) is completely specified 
by its second moment, it is sufficient to describe the corre- 
lation function (Tr AH x TiBTL) for two arbitrary sources 
A and B. The correlation law we choose is 

(TyATL x TrBTL)/v 2 
= Tr(A - E X A T E X ) \b+(1- e t )rB T T- x 



+(l-e s )Y,JkBJ k (8) 

fe 

-(l-e.)(l-et)X) JkTB T T~ 1 J k . 



For e s = e t = the correlation law is invariant under both 
reversal of time and rotation of spin. This is the symme- 
try class CI. A nonzero value of breaks time-reversal 
symmetry. Therefore, by increasing we cross over to 
class C. A nonzero value of e s breaks spin-rotation in- 
variance, so by increasing e s we cross over to class Dill. 
By increasing both e s and et we break all symmetries and 
cross over to class D. We call a symmetry "maximally 
broken" when its symmetry-breaking parameter (e s or 
et) equals unity. Whenever a symmetry is either unbro- 
ken or maximally broken, the probability distribution of 
the Gaussian ensemble can be presented in the simple 
form (Q), with the corresponding symmetry constraints 
imposed on TL. 

All information about the level statistics is contained 
in the joint probability distribution for the eigenvalues, 
P({uj}). This distribution is a complicated function in 
general, but it takes a simple form for each universality 
class. By diagonalizing the BdG-Hamiltonian 



TL = U 



uj 
-uj 



diag(cJi,w 2 , ...) 



and computing the Jacobian of the transformation to di- 
agonal form, we obtain the formula 

P({Lu})d{uj} - [J \lo 2 - u 2 f ft \u k \ a e-"ll* 2 du> k (9) 

i<j k 

where, for the individual cases, 

class D : (3 = 2, a — 0; 
class C : (3 = 2, a = 2; 
class Z3III : (3 = 4, a = 1; 
class CI : (3—1, a = 1. 

These expressions for P({uj}) can be derived by elemen- 
tary means. A particularly elegant derivation uses the 
interpretation of TL as being tangent to the symmetric 
space of type D, Dill, C, or CI. The Jacobian can then 
be read off immediately from the tabulated root systems 
of these spaces. 

The formula for P({oj}) permits some immediate con- 
clusions to be drawn. Clearly, the significance of the pa- 
rameter a is that it governs the level repulsion from the 
origin uj — 0, while (3 gives the mutual repulsion between 
levels. For the following it is useful to view the factor 
\ujk\ a as being due to the interaction of the k th level with 
its "image" at —uJk- Similarly we view the factor uji + ujj 
in ujf — uj 2 as resulting from the interaction of the i th 
level with the image of the j th level. At energies uj much 
greater than the mean level spacing, the interaction of 
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levels with their distant images at negative energies is 
expected to be irrelevant. Therefore the level statistics 
derived from (||) will reduce, in that limit, to the usual 
Wigner-Dyson statistics as determined by the parameter 
[3. On the other hand, in the opposite limit of energies 
of order unity on the scale set by the level spacing, the 
level statistics will be different from Wigner-Dyson. In 
particular, by the definition of P({w}) as a joint prob- 
ability density we immediately conclude that the mean 
density of levels near zero behaves as 



<pH> = (Tr8(co - H)) ~ M° (w -> 0). 



(10) 



Note that for the systems where our random-matrix de- 
scription applies, the exponent a is predicted to be uni- 
versal, dependent only on symmetry. The value of a for 
the symmetry classes CI and C is easily understood from 
the fact that the repulsion of a level from its own image is 
caused by the pairing field A. For class C pairing matrix 
elements are complex, whereas for class CI all pairing 
matrix elements can be chosen to be real. By a stan- 
dard power counting argument this results in a — 2 and 
a = 1 respectively. To understand why a is zero for class 
D, note that in this case a level and its own image are 
not really physically distinct but are copies of the same 
single-electron state. (In contrast, for the classes C and 
CI the hole level has its spin flipped relative to the parti- 
cle level.) The pairing matrix element between identical 
states vanishes by the Pauli principle - or put differently, 
the matrix A in (||) is antisymmetric and therefore has 
zeroes on its diagonal -, which results in a = 0. 



V. SPECTRAL STATISTICS 

A. Exact Results 

Our interest here is in the level correlations for a large 
matrix dimension 2N. These are easy to compute when 
the symmetry class is C or D. Consider first class C. 
For this symmetry class we can interpret P({lo}) as the 
joint probability density of a Gaussian Unitary Ensem- 
ble (GUE) of 2N levels w\,u}\, ...,ujn,^n subject to the 
mirror constraint lo^ = —cut- The GUE joint probabil- 
ity density, in turn, can be interpreted as the square of 
the ground state wavefunction for a system of spinless 
nonrelativistic noninteracting ID fermions confined by a 
harmonic well jlffl. This correspondence of levels and 
Fermi particles turns the n-level correlation functions of 
the GUE into the n-point static density correlation func- 
tions of the Fermi system. In the large- N limit, the spa- 
tial variation of the harmonic confining potential becomes 
(locally) negligible and the gas of fermions can be treated 
as locally free. The mirror constraint means that when- 
ever a fermion approaches zero, then so does its mirror 
image. Because the Pauli principle makes the wavefunc- 
tion vanish as two fermions approach each other, this 
amounts to hard wall (or Dirichlet) boundary conditions 



at to = 0. Hence we can compute the level density and 
its correlations as the particle density and its correla- 
tions for a free ID Fermi gas with Dirichlet boundary 
conditions at the origin. The free-fermion wavefunctions 
that vanish at lu — are sin(cJr), where r plays the role 
of a "wavenumber" . By summing over the Fermi sea of 
states occupied in the ground state, we obtain for the 
mean density of levels 



sin 2 (wr)dT = 1 



sin 2ttlo 
2ttlo 



(11) 



Here and throughout this subsection we follow the con- 
vention of measuring u> in units of the mean spacing be- 
tween neighboring particles (i.e. of the level spacing) at 
a distance of many spacings from zero. Note that (p(oj)} 
for lu — > has the behavior expected from ( |l0| ) (recall 
a = 2 for class C). A similar calculation of the density- 
density correlator of the Fermi gas yields the two-level 
cluster function: 



(p(wi)p(w 2 )) - (p(wi))(p(w 2 )) 
- [8{oji - lo 2 ) + 5(u>i + lo 2 )] (p(wi)) 
sin7r(wi — L02) sin7r(aji + L02) 



Tt(lVi — LO2) 



7r(cJl + LO2) 



Keeping r = L0\ — 102 7^ fixed and letting u>\ + lu 2 tend 
to infinity, we recover the familiar GUE two-level cluster 
function — sin 2 (7rr)/(7rr) 2 . Similarly, all n-level functions 
i? n (wi, to n ) = (p(u!i)...p(Lo n )) can be calculated. On 
subtracting the level self-correlations, which amounts to 
normal ordering in the particle-gas formulation, we ob- 
tain the result 



R n (tu 1 ,...,w n ) =Det[*c(wi,Wj)]j,j=i,..., n , 
2 f* 

^ ci^ii^j) = — / sm(ujiT) sm(ujjT)dT, 
* Jo 



(12) 



by simply using Wick's theorem for the free Fermi gas. 

We turn to the symmetry class D. It is convenient 
again to use the interpretation of the joint probability 
density as a Gaussian Unitary Ensemble of 2N levels 
with a mirror constraint. The only change from before is 
that the repulsion of a level from its own mirror image is 
now absent (a = 0). Correspondingly the single- fermion 
wavefunctions of the Fermi gas no longer vanish on ap- 
proaching the origin. Instead, what we need to demand 
is that they be even functions of to, which is the same 
as imposing vanishing derivative (or Neumann) bound- 
ary conditions at to = 0. Thus the level «-> particle corre- 
spondence now leads to the free Fermi gas with Neumann 
boundary conditions at the origin. Doing the same kind 
of calculation as before we find 

, / m 2 f v o, „ , sin(27uj) . , 
(p(u>)) = - cos 2 ( W r)dr = l + ^ i, (13) 

and the result ( |l2| ) remains valid if we replace by 'I'd, 







(wi , ujj ) = — / cos(wir) cos(cJjr)rfr. 



From (13) we see that for a metallic quantum dot with 
spin-orbit scattering (class D), the proximity of a super- 
conductor with {e 1 ^) — enhances the level density at the 
chemical potential! While this effect may seem physically 
surprising, it is very natural in the Fermi-gas picture of 
the levels. The pressure of the gas pushes particles (or 
levels) against the "wall" at ui = 0. Because it is the 
current rather than the density that is required to vanish 
by the Neumann boundary condition, an excess particle 
density forms at the wall such that the extra statistical 
force balances the pressure. 

More effort is required by the symmetry classes CI 
and Dili, where (3=1 and (3 = 4. It is still possible 
in these cases to map the level statistics problem on a 
model of particles moving on a half line, but progress is 
slowed down by the fact that the particles now interact 
with each other. By a standard transformation Jl8[ | one 
can show that their motion is governed by the Hamilto- 
nian of the Calogero-Sutherland model (CSM) associated 
with the symmetric spaces of type CI and DHL For the 
case of the CSMs corresponding to the Wigner-Dyson 
Ensembles, it was recently found |n| that the CSM par- 
ticles behave as a gas of free anyons, i.e. particles with 
fractional charge and statistics. Although we have some 
preliminary results indicating that the free anyon gas pic- 
ture can be adapted to the present situation, the details 
have not been worked out yet. 

A quick way to get the infrared (or large-w) asymp- 
totics of the level density for CI and Dili is to bosonizc 
p0| the CSM. This procedure has been argued Q to 
yield the c = 1 conformal field th eory of a free boson 
with compactification radius R = \J (3/2. The expression 
for the CSM particle density rpip in terms of the boson 
field <p is H 



i/np = d u (p + const x cos(v47r<y9/i? + kpu). (14) 

(Recall that by the level «-> particle correspondence w 
is to be interpreted as a space coordinate here.) The 
mirror constraint of the CSM for CI and Dili translates 
into a boundary condition on tp at lj = 0. Since the ver- 
tex operator exp(>/4ni<p / R) has the scaling dimension 
1/R 2 = 2/(3, we expect 



( P ( w )) = i + w - 2 /%H + 



(15) 



where Ap(u>) is a function that oscillates with a period 
determined by the mean spacing. Note that (15) is con- 
sistent with the (3 = 2 results Jll] ) and (|l3|). Note also 
that the first term on the right-hand side of (14) gives a 
vanishing contribution to the average density, although 
it does contribute to the CSM density-density correlator. 
This is because the current d u ip, being linear in the bo- 
son field ip, has a vanishing expectation value even when 
there is a boundary. 



The validity of bosonization and conformal field theory 
arguments is restricted to the infrared regime. To obtain 
expressions that are valid in the entire range of frequen- 
cies io, we turn to the orthogonal polynomial method of 
Dyson and Mehta ^3|. The substitution Xk = w| turns 
(y) for a = 1 into 

p({x})d{x} = const X T\\xi - XjfT\e~ Xh / v dxk, 

i<j k 

which defines what has been called |24[] the Laguerre Or- 
thogonal Ensemble (LOE) for (3 = 1, and the Laguerre 
Symplectic Ensemble (LSE) for (3 = 4. Note that this 
nomenclature is rather unfortunate in the present con- 
text. As we saw, the LOE relates to the symmetric 
space Sp(2iV)/U(iV), while the LSE relates to the sym- 
metric space SO(4JV)/U(2AT). In both cases the invari- 
ance group is a unitary group, U(A) or U(2AT). Closed 
expressions for the n-level correlation functions of these 
ensembles have recently been published by Nagao and 
Slevin |23] . By using some identities for Bessel functions 
and returning to the variable 10 = x 2 , we can cast their 
results for the mean density in the form 



CI 
Dili 



(p(w)) = F(ncj), 

(p(co)) = F(2ttuj) + 7tJi(2tuj)/2, 

7T 



F(z) 



dt J (i)Ji (*)/*, 



where Jk is the Bessel function of order k. (Remember 
that we are taking the level spacing at large u> for our en- 
ergy unit. The levels are counted without multiplicity.) 
From this we read off the small-w expansions: 

(p(co)) = (3n 2 uo/A + O{uo z ) (CI and Dili). 

Knowing the mean density, we can construct the full one- 
point function (g(to)) = (Tt(lu + iS — H) _1 ) by causality, 
i.e. by using the dispersion relation that connects the real 
and imaginary parts of a holomorphic function on the up- 
per complex half-plane. The results can be presented in 
the form 

CI : (g(co)) 

/oo f+1 _ y 2 gi7rw(«— v) 

du / dv 



Vu 2 - 1 u—v 



Dill : (g{w)) = —in + in du y 

Ji \Ju 2 — 1 



a 



217TLJU 



' , f +1 , \At 2 - 1 e 2inu( - u -^ 

-i I du / dv —== . 

/i J-i \Jl-v 2 u-v 

Although it is hard work to construct these expressions 
directly, they can easily be verified. For that we sim- 
ply differentiate the result for (g(u>)) with respect to u>, 
thereby cancelling the factor u — v in the denominator 
of the double integrals. The integrals over u and v then 
separate, and on taking the imaginary part each integral 
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produces a Bessel function. By using standard recur- 
sion relations for these functions and then undoing the 
(^-differentiation by integration, we immediately retrieve 
the expressions for (p(ui)} given earlier. 

The double integrals for CI {(3 = 1) and Dili [fi = 4) 
are seen to be related by a duality transformation that 
exchanges the compact (v) and noncompact (u) degrees 
of freedom. A similar duality relation holds for the con- 
ventional Wigner-Dyson ensembles with (3—1 and (3 = 4 
H . In the limit of large u> we get the following asymptotic 
expansions for the one-point function: 



e 



Dill : (g(co)) = -in + ±- + ™e 2 ™ + ™/ 4 + 

ALU 2y/UJ 

For completeness, the one-point functions for the sym- 
metry classes C and D (/3 = 2), as determined from ([l - 
and dl^ ) by causality, are 



(sM) = - i7r + (i - a )- 



J2-KIUJ 



2uj 



(C and D). 



By comparing with ( |l5| ) we see that the oscillatory cor- 
rection to the stationary asymptotic limit (g(u>)) — * — itr 
agrees with what is expected from the conformal limit of 
the Calogero-Sutherland model, in all cases. The smooth 
(1/cij) part of the correction is purely real and does not 
enter into the asymptotic expression for the density of 
states. 

The authors of Ref. 0| subjected the n-level cor- 
relation functions for n > 1 to a renormalization or 
unfolding procedure in the low-frequency regime they 
call "nonuniversal" (meaning different from standard 
Wigner-Dyson). We wish to emphasize that such a pro- 
cedure is neither necessary nor appropriate here. Both 
the mean density and the level correlation functions are 
universal as they stand - the restriction to the class of 
system we have delineated being understood, of course - 
and are not to be corrupted by any kind of unfolding. 



B. Diagrammatic Perturbation Theory for the 
One-Point Function 



As we have seen, (g(ui)) tends to a constant for fre- 
quencies much larger than the level spacing. The leading 
smooth (i.e. nonoscillatory) correction is of order 1/uj 
in all cases. More precisely, on restoring the physical 
units and taking into account the multiplicity of levels, 
we have (g(w)) = —iirv + c/w + O(l/oj 2 ), where c = —1 
for C and CI, and c = +1/2 for D and DHL v is the 
asymptotic (i.e. large-w) density of states. We are now 
going to show how to obtain this result by a variant of 
the impurity diagram technique, a method which has the 
attractive feature of generalizing easily to the calculation 
of transport properties of an open system. It also has the 



great virtue of lending itself to semiclassical interpreta- 
tion, which will help improve our understanding of the 
physics involved. 

The impurity diagram technique in its present version 
starts from the usual idea of expanding (uj + iS — 7i) ~ 1 in 
a geometric series with respect to Ti and then taking the 
ensemble average. Because TL is Gaussian distributed, 
the ensemble average is evaluated by forming all products 
of pairwise contractions (TCH) , which are determined by 
the basic law (|J). To resum the relevant contributions, 
we use standard diagrammatic techniques. On multiply- 
ing the factors on the right-hand side of (||) we generate 
eight terms. In explicit index notation these are given by 



nf)„ s =s a5 6. 



(IT?) 



0:7, /36 



(n c i) Q 
(n c .°) Q 



■y,/38 



Tq 7 (t )spi 
y^(JfcT)a 7 ( 



_1 >/fe)(5^, 



= — / y (Jk^x)jaC^x Jk)f35, 



(16) 
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It is characteristic of the contractions indexed by the let- 
ter "c" that the initial states (3, S bear a definite relation 
to each other, and so do the final states a, 7. This sit- 
uation is reminiscent of the cooperon channel of disor- 
dered mesoscopic systems where a pair of particles with 
initial momenta p and —p are scattered to final mo- 
menta p' and —p'. Similarly, the contractions indexed 
by "d" correspond to the diffuson channel where a pair 
with momenta p, p' is scattered to a pair with momenta 
p' \p. The contractions with subscript "D" owe their 
existence to the operation of particle-hole conjugation 
X 1 > —Tj x X t Yj Xi whose fixed point set is the orthogonal 
algebra D2N = so(4iV). The name of the "^4" -type con- 
tractions is motivated by the fact that they determine 
the second moments of the conventional Wigner-Dyson 
ensembles describing N-systems (without any coupling of 
particles and holes), which derive from the unitary alge- 
bra Ak-i = su(fc). The numerals and 1 distinguish 
between spin-singlet and spin-triplet contractions. Using 
the conventions ( jig ) we can write the correlation law (@) 
in the form 



n= n 



dO 



n*; 



-(i-e t )(ni° + n^) 

+ (l-e s )(l-e t )(U c 2+K 1 ) 



(17) 
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Our goal is to find the large-w behavior of (g(u>)). What 
are the dominant diagrams in this limit? From what has 
been said, the A-type contractions give rise to Wigner- 
Dyson statistics, whereas the D-type contractions are re- 
sponsible for the corrections to it. Since our systems are 
Wigner-Dyson in the limit of large to, the D-type contrac- 
tions must become irrelevant in that limit. Moreover, 
in the Wigner-Dyson regime the average Green's func- 
tion is known to be featureless and independent of the 
symmetry-breaking perturbations e s and e t . We there- 
fore conclude that (g(to)} is completely determined by 
IT^ -contractions for lo — > oo (and large N). By sum- 
ming all nested self-energy graphs, we get Pastur's 
approximation G° to G = ((u + iS — TL)^ 1 ): 



G° = (lj + id - v 2 TrG°Y 



(18) 



This equation is exact for N — > oo and large lo. Its so- 
lution yields Wigner's semicircle law for the density of 
states. Putting v 2 — X 2 /AN and focusing on the central 
region of the semicircle, we obtain 

TrG° = — vkv + {ttv) 2 lo/8N + 0(lo 2 /N 2 ) 

where v = AN/ttX is identified as the asymptotic density 
of states. What we need to do to probe the local structure 
of the spectrum, is to keep the product vlo fixed while 
sending N to infinity. The corrections to TrG° = —iirv 
from Pastur's equation are seen to become negligible in 
this limit. However, we know that corrections to the 
stationary asymptotic behavior TrG° = —iirv do appear 
as we approach zero frequency. These must be due to 
the contractions of type D. The leading correction is 
depicted in Fig. ^, where the light-shaded regions repre- 
sent ladder diagrams built either from ITp-contractions 
or from IT^-contractions. 



+ 



FIG. 2. Diagrams contributing to the average sin- 
gle-particle Green's function. The light-shaded regions rep- 
resent a D-type cooperon mode, the dark-shaded one a non- 
singular n^°-ladder. 

The sum of the former diagrams, which we call the re- 
type spin-singlet cooperon mode and denote by 5° p S , 
has the expression 



5° = v 2 U /(l - v 2 KUq) 



with K ai ^ s = S af3 S 7S G° aa G ^ and n = 11$. Its large-iV 
limit is 



5° = 



A 2 n 



-ITTVLO 



<D(l/N). 



(19) 



Similarly, the sum of all ladder diagrams built from 11^- 
contractions, the D-type spin-triplet cooperon, is evalu- 
ated as 



s 1 = (1 - e s )^ 2 n x /[l - (1 - eJ^MJ 

A 2 n x 



Vs 



Oil /N) 



(20) 



where rj s — ANe s and ID = 11^}. The dependence of 
S 1 on the parameter e s through the product iVe s means 
that the breaking of spin-rotation invariance takes place 
on scales e s ~ 1/N and thus is very fast. This rapid 
crossover happens because the crossover scale is deter- 
mined by the typical size of a symmetry-breaking ma- 
trix element in relation to the level spacing, which is 
v~ x = Xir/AN for our choice of normalization. 



Note that the expressions ([L9|) and ( |20D are singular 
at lo = = r/ s even though the sums of ladder diagrams 
they represent are built from retarded Green's functions 
only (G + G + channel). This is a novel feature which 
does not occur for the standard Wigner-Dyson ensem- 
bles, where singular ladders exist only in the advanced- 
retarded (or G~G + ) channel. The singularity in the 
present case comes about because the minus sign from 
K ai ^ ai = —1/X 2 is cancelled by a minus sign appearing 
in the definition of the contractions of type D, thereby 
turning an alternating (conditionally convergent) series 
into a divergent one. 

The dark-shaded region appearing in the second di- 
agram of Fig. U represents a Il^-ladder. According to 
( jig ) , the contractions of type A come with an overall plus 
sign, so the minus signs now do not cancel, and the ladder 
sum is always finite. Computing the sum we find that this 
nonsingular IT^-ladder renormalizes the first diagram in 
Fig. | by a factor of 1 - (1 - 1 + 1 - ...) = 1/(1 + 1) = 1/2. 
(We mention in passing that nonsingular ladders of this 
kind are the random-matrix analog of the single- impurity 
lines that appear in the context of the impurity diagram 
technique.) 

To evaluate the first diagram of Fig. || we need the 
following sums: 



E^V^-E^W^ 



'/3a 





E (ns: 



a/3, /3a 



E ( S ^)/3a (Jk^x)p a - +3. 





We then obtain 

(g(u>)) = — iTTV + %-kv 



3/2 



1/2 



77s — TTIVLO 



which agrees with the result stated at the beginning of 
the current subsection for rj s — (classes G and GI) and 
■q s — » oo (classes D and Dili). For the classes G and D, 
smooth corrections of higher order in 1 /lo are completely 
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absent from the exact result of Sec. VA. This implies 



that all diagrams of higher order than the ones consid- 
ered here must cancel each other, in these two cases. (The 
oscillatory correction ~ Apuj^ 2 / 13 is nonanalytic in the 
expansion parameter 1 /lu and therefore remains invisible 
to all orders of perturbation theory.) 



VI. SLOW MODES 



In all previous Green's function treatments of NS- 
systems the diagrams were enumerated by the number of 
Andreev reflections. Unfortunately, when the perturba- 
tion expansion is organized in that way, the vast number 
of possibilities to insert Andreev reflections into the dia- 
grams generates a flood of terms which is hard to control, 
and as a result it is very easy to miss important contri- 
butions. The technical innovation made in the present 
paper is not to single out Andreev reflections but treat 
them on exactly the same footing as the processes of im- 
purity scattering. This is possible by our dynamical as- 
sumptions ensuring that the quantum mechanical phase 
acquired during Andreev reflection, can be regarded as a 
random variable with zero mean. Our key technical step 
is the decomposition ( |l7| ) which leads to an organization 
of the perturbation-theory diagrams by symmetry. In the 
preceding subsection we discussed how the contractions 
H c £ and generate singular geometric series of ladder 
diagrams. In the same way, every one of the other con- 
tractions gives rise to one singular ladder. These singular 
modes can be visualized as follows: 



P(h) 
P(h) 



G+(G" 



P(h) 
P(h) 



G~(G+) 
A— type cooperon, 

G+(G~) 



P(h) 
h(p) 



G-(G+) 
A— type diffuson, 

G+{G-) 



G+{G-) 
D— type cooperon, 



P(h) 
h(p) 



G+(G~ 



G+{G~) 
D— type diffuson. 



The dotted vertical lines represent both impurity scat- 
terings and Andreev reflections, and they denote any 
one of the eight contractions IF^ (X = A, D; x = c, d; 
S = 0,1). The type of contraction is invariant within 
one ladder. The A-type modes are built from states of 
identical charge (two BdG-particles or two BdG-holes) 
propagating on opposite segments of the ladder, whereas 
the D-type modes are built from charge-reversed states 
(one particle and one hole). The former are singular in 
the G+G- channel, the latter in the G+G+ (or G~G~) 
channel. The arrows on the Green's function lines indi- 
cate the order in which single-particle states are visited. 
For the cooperon modes the order on both lines is the 
same, while for the diffuson modes it is reversed. In the 
limit u> = ?7 S = rjt — (with rj t = 4iVet) all modes are sin- 
gular, or massless. The -D-type modes are made massive 
by frequency (or voltage) u> while the A-type modes are 
insensitive to such a perturbation. The A-type cooperon 
and the D-type diffuson are made massive by the break- 
ing of time-reversal symmetry. Since a Green's function 
line carries spin 1/2, the modes decompose into spin- 
singlet and spin-triplet ones. The spin-triplet modes are 
sensitive to spin-orbit scattering while the spin-singlet 
modes are not. 

We wish to mention that there is some redundancy 
in our classification of modes, as the basic particle-hole 
symmetry (^) causes the existence of certain relations 
among the matrix elements of the Gorkov Green's func- 
tion G ± (w) = (uj±ie — TL)^ 1 . In particular, the particle- 
particle and hole-hole matrix elements are related by 



G±M = -GU-") T - 



(21) 



Similarly, Gp h (w) = — G^ p (— lu) t . These identities tran- 
scribe into relations connecting the singular modes. For 
example, by using ( pl| ) on one of the Green's function 
lines of the D-type cooperon, we can make this mode 
look like the A-type diffuson, at the expense of having 
to change the sign of one frequency (lu —> — lu). In a 
similar way, the D-type diffuson is related to the ^4-typc 
cooperon (again with a sign change in one of the fre- 
quencies). In spite of that, we prefer to treat the A- and 
D-type modes as separate entities. The main reason for 
doing so is that they respond differently to translations 
of the energy: while the D-type modes are made massive 
by shifting the energy, the A-type modes are not. 

In the present paper we restrict our considerations to 
the ergodic (or zero-dimensional) limit. To go beyond, 
we should associate with each A-type and D-type mode 
a small momentum variable ( "slow modes" ) and sum over 
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momenta. In this way, it will not be difficult to generalize 
our results beyond the ergodic limit. 



VII. WEAK LOCALIZATION 

Having made a thorough analysis of the isolated An- 
dreev quantum dot, we now turn to the discussion of the 
associated open system and its transport properties. To 
open up the dot in the simplest possible way, we couple 
it to a single lead with 2M open channels (the factor of 
2 accounts for the spin degree of freedom), see Fig. 0. 




FIG. 3. Andreev quantum dot (N) coupled to a single lead 
(L) via a tunnelling barrier (T). The flux loop on the right 
is introduced to adjust the difference of the order parame- 
ter phases of the superconducting re gion s (S) to the value 
(j>i — <f>2 = 7T (cf. the discussion of Sec. III). 



The transmission of charge excitations from the lead 
to the interior of the dot is modelled by a set of (spin- 
independent) real hopping matrix elements W^ a , where 
the index a = 1,...,M (/i = 1, . . . , N) enumerates the 
channels carried by the lead (the sites of the dot). We 
assume AT > M > 1. Let g denote the conductance 
measured in units of e 2 jh. To calculate g, we employ the 
Landauer-Buttiker formula as generalized to NS-systems 
by Takane and Ebisawa |25|, 



!1 



hpi2 
la 1 ' 



(22) 



where the composite label a — (a, s a ) comprises the spin 
s a = ±1/2 and the index a of an open channel in the 
lead, and S~? denotes the scattering amplitude connect- 
ing a particle coming in in channel a with a hole going 
out in channel b. The S'-matrix is given by 



where 



G=(iS-H + iWW T ) 



(23) 



is the Gorkov Green's function evaluated at the chemical 
potential. Without loss of generality, we may assume the 
matrices W — {W^a} to be of the form 



W^ a = 7 1 / 2 (5 AIQ ( M = 1, ■ • • , N; a = 1, . . . , M). (24) 

The unitary transformation necessary to transform W to 
the form ( p4f ) can be absorbed in the Hamiltonian H by 
the invariance properties of the random-matrix ensemble 
(|). Combining Eqs. ( p3| ) and ( p2| ) and making use of 
( 24|) we obtain 

3 = 8 £ r M |G(^ Sih ),( M /v iP )| 2 r^, (25) 



where 



7 for fi < M, 
else. 



We are going to calculate this expression to leading order 
in the small parameters l/N, M/N and next-to-leading 
order in 1/M ^6|. Owing to the presence of the BdG 
particle-hole degree of freedom, an analysis of Eq. fl25| ) 
within the framework of plain diagrammatic perturba- 
tion theory turns out to involve a sizable number of dia- 
grams. It is more efficient to pre-analyze j25|) by means of 
a set of exact identities (Ward identities) before turning 
to diagrammatic methods. This calculation is detailed 
in Appendix Here we restrict ourselves to a presen- 
tation of the results and their interpretation in terms of 
semiclassical trajectories. 

Schematically, the conductance can be represented as: 
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where the wavy lines stand for the quantities {F^}, the 
shaded region denotes the singular A-typc diffuson mode 
introduced in Sec. [v| , 



and a summation over indices is understood. The WL- 
building block represents a quantum interference correc- 
tion (the NS-analog of the well-known weak localization 
correction for normal metals) to the classical conduc- 
tance. In contrast with the pure N-case, however, two 
qualitatively different processes contribute to the weak 
localization correction for the Andreev dot: 



WL 



A + 



D 



Here the A(D)-block is due to the presence of singular 
A(Z))-type modes. Whereas the A-type contribution re- 
sembles the standard weak localization correction known 
from normal metals, the D-term does not have any ana- 
log in pure N-systems and is of a different nature. In the 
following we discuss separately the classical conductance 
(the first diagram in (|2li|)), the >l-type correction, and 
the D-type correction. 
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Classical conductance: Qualitatively speaking, the 
conductance is given by 



(27) 



where A^ is the amplitude to traverse a certain scattering 
sequence (indexed by i) connecting an incoming particle 
channel with an outgoing hole channel. The classical 
value of the conductance, go, is obtained by evaluating 
the incoherent sum 

i 

Quantitatively, we obtain 

5o = 2MT, 

where the transmission coefficient T is the probability 
for an electron incident from the lead to enter the dot 
instead of being reflected back into the lead ]27| ]. This 
result is easy to understand. By the ergodicity of our 
system, an electron leaves the dot with equal (1/2) prob- 
ability as a particle or as a hole. In the latter case, two 
elementary charges are transferred across the entire sys- 
tem. Thus the dimensionless conductance per channel is 
2 x 1/2 x T = T. Multiplying by the number of channels 
we get go = 2MT. 

A-type corrections: Weak localization corrections to 
the classical conductance originate from the phase- 
coherent contributions of nonidentical paths to the sum 
of amplitudes ( |27| ) . In the case of the A-type correction, 
such contributions are due to pairs of paths that differ 
by a sequence of scattering events traversed in opposite 
directions as is indicated in Fig. [|. 




FIG. 4. Pair of semiclassical paths contributing to the 
A-type weak localization process. 



The sum of these "maximally crossed" segments of pairs 
of paths is represented by the building block A in (26). 
More specifically, 



shown merely for the sake of clarity but do not contribute 
to the A-block as such, and the dots stand for diagrams 
of a more complex structure that have to be taken into 
account to obtain a result consistent with unitarity. It is 
the presence of these unitarity-preserving contributions 
that renders the calculation of the A-block within plain 
diagrammatic perturbation theory lengthy. The alter- 
native computational scheme presented in Appendix ^ 
yields 

y 2 \MT + r) t MT + Va + T]t 

where the parameters n s = 47Ve s and rjt = 4iVet are the 
scaled symmetry-breaking parameters of our model. 

D-type corrections: A pair of paths contributing to the 
-D-type weak localization process is shown in Fig. ^. 




FIG. 5. Pair of semiclassical paths contributing to the 
D-type weak localization process. The triangles represent An- 
dreev reflections. 

Note that the self-intersecting loop must contain a non- 
vanishing even number of Andreev reflections (the figure 
displays the simplest possible case of just two Andreev 
events). We note in passing that the A-type loop shown 
in Fig. ^ may contain Andreev reflections, too (for this 
reason we said that the NS A-type correction is analogous 
to, though not identical with, the normal weak localiza- 
tion correction), their presence is just not imperative like 
in the D-case. Clearly, the D-type correction does not 
have any analog in normal metals. Note also that the 
closed loop in Fig. || involves only one of the two paths. 
This shows that the existence of the D-type process is 
essentially due to the nontrivial behavior of the single- 
particle Green's function. The same mechanism of quan- 
tum coherence at the single-particle level was responsible 
for the correction to the single-particle density of states 



discussed in Sec. [VB . 

In diagrammatic language, the loop insertion in Fig. |^ 
is represented by 



D 



where the shaded region represents an A-type cooperon, 
the subscript 't' means that the external arrows are 



t 

where the shaded region now represents a D-type 
cooperon mode and the dots stand for a set of unitarity- 
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preserving counter diagrams. The quantitative analysis 
yields 



Sg D = (1 - T) 1 - 



3MT 

mtTvs 



A striking feature of this expression is its insensitivity 
to the breaking of time-reversal symmetry: the D-type 
weak localization correction for NS-systems survives the 
application of external magnetic fields f28|| . Collecting 
terms we obtain the final result 



(g)=2MT- 



1 - T) ( 1 
1 



M 2 / 
T 

2 + % 



3MT \ 
MT + rjs) 

3 ) 

+0(l/M,M/N) 



(28) 



for the dimensionless mean conductance of our system at 
zero bias. We see that the D-type correction is zero for 
T = 1, while the A- type correction vanishes in the limit 
T — > 0, with MT held fixed. From what has been said 
about the D-type modes, we expect the D-type correc- 
tion to disappear at finite bias. Detailed analysis shows 
that the crossover scale for this to happen is determined 
by vuj ~ MT. 



VIII. UNIVERSAL CONDUCTANCE 
FLUCTUATIONS 

The conductance fluctuations of normal-conducting 
systems |2j| have been studied extensively. They are 
independent of system size and strength of the disorder 
and depend only on symmetry. The latter dependence 
can be summarized by saying that var(g) is proportional 
to the number of massless modes for a given universal- 
ity class. When all symmetries are broken, the ^4-type 
spin-singlet diffuson is the only mode which is massless. 
As we switch off the spin-orbit interaction, the A-type 
spin-triplet diffuson modes become massless, too, which 
increases var(g) by a factor of four. If in addition time 
reversal is a good symmetry, the cooperon modes become 
massless, thereby increasing var(g) by yet another factor 
of two. 

The NS-systems considered in the present paper are 
expected M to show conductance fluctuations that are 
qualitatively similar to those of N-systems. To calculate 
the variance, we may use an extension of the diagram- 
matic method described in the previous section or, al- 
ternatively, we may map our random-matrix model on 
a zero-dimensional field theory of the nonlinear a model 
type. In the present paper neither of these methods will 
be used. Instead, we will turn to another approach, which 
is restricted to the strong-coupling limit T = 1 but has 
the great advantage of being very simple. 

The symmetry properties of the S-matrix derive from 
the symmetries of the Hamiltonian by exponentiation. 



As before, let M denote the number of channels in the 
lead, not counting spin and particle-hole degeneracy. By 
the considerations of Sec. [n] the S'-matrix may be re- 
garded as an element of the symmetric space SO(4M) 
for class D, Sp(2Af) for C, SO(4Af)/U(2Af) for Dili, 
and Sp(2M)/U(M) for CI. We will refer to these spaces 
as "S'-matrix manifolds" for short. Let A = (a, s, a) 
(a = l,...,M;s — ±1/2; a = p,h) be a composite 
index. From the definition of the transmission coeffi- 
cient T as a "sticking probability" |u| we have T — 
1 - \(Saa)\ 2 + 0(1/MT). Therefore T = 1 implies an 
S-matrix with vanishing ensemble average, which means 
that S can be taken to be uniformly distributed on its 
S-matrix manifold. 



A. Class C 

For the symmetry classes C and CI the S-matrix oper- 
ates on the tensor product of channel space and particle- 
hole space, while spin is accounted for by multiplication 
of the conductance by a factor of two. Recall the defini- 
tion of the symplectic group Sp(2M) by 



u~ 1 ^ = u = cu~ lT c- 1 



(29) 



where C = 1m <8 i<J y . In keeping with the above, we 
take the S-matrix S = U for class C to be uniformly 
distributed on Sp(2Af). In other words, ensemble aver- 
ages (...) are computed by integrating with respect to the 
Haar measure dU: 



(f(U)) = 



Sp(2M) 



f(U)dU. 



The canonical projection of Sp(2M) onto the coset space 
Sp(2M)/U(M) by U i * UU T turns the Haar measure of 
the former into the invariant (or uniform) measure of the 
latter. Therefore, ensemble averages for class CI can be 
obtained from 

(f(S)) C i - (f(UU T )). 

Because the Haar integral is invariant under left and right 
translations, 



f(U L UU R )dU 



'Sp(2M) JSp(2M) 

the defining equations for Sp(2M) lead to 

(Uab) = 0, 
{U ab U*cd) 



f(U)dU, 



(UabUcd) 



S AC S BD /2M, 
C AC Cbd/2M. 



(30) 



To compute the ensemble average of a product of two 
U's and two £/*'s we note that if ipA & V are the com- 
ponents of a vector transforming according to the funda- 
mental representation of Sp(2M), there exist only two 
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independent invariants on V <g> V* <E> V ® V*, namely 

J^^a^a^b^b and Y, c abCcd^a^b^* c ^*d- Using this 
elementary group-theoretical fact we obtain 



(^AlBi ^CiDi ^A 2 B 2 ^C 2 -D 2 

2M- 1 



2M(2M+1)(2M- 



1 

'2M(2M+ 1)(2M- 2) 



^d^CjfetD^B, D 2 

+<5a 1 c 2 ^ 2 Ci<5b 1 d 2 ^b 2 _Di 
+C^ lJ 4 2 Cc lC2 C SlS2 C £ ) 1 £) 2 

^iCi^aCa^BiUa^Oi (31) 

-<5aiC 2 ^a 2 Ci <5BiDife 2 n 2 



+(<5a 1 Ci^a 2 c 2 - ^ 1 c 2 ^ 2 c 1 )Cb 1 b 2 Cd 1 l> 2 

+CAiA 2 CciC 2 (<5BiDife 2 D 2 — 8b 1 D 2 ^B 2 D 1 , 



The numerical coefficients in this expression are deter- 
mined by summing over any two pairs of equal indices 
and then comparing the results to ( |30| ) using the rela- 
tions (p9|). Equation ( |3l| ) entails 

E^p^^b^b) = (2M + l)" 1 , 

AB 

which can be used to compute the weak localization cor- 
rection for class CI. Summing over initial and final (or 
particle and hole) channels, we get (TrS' ph S' thp ) C i = 
M 2 /(2M + 1) which yields Sg — —1 in agreement with 
(p8|). To calculate the conductance fluctuations for class 
C, we deduce from risi]) 



((Tr S ph S thp f 



Subtracting the square of the first moment and multi- 
plying by a factor of 4 x 4 for charge and spin, we get 
var(g) = 2. 



B. Class D 

The symmetry class D can be treated by direct tran- 
scription from class C, the only difference being the way 
the spin enters. The ^-matrix now operates on the full 
tensor product of channel space, particle-hole space and 
spin space. The S'-matrix manifold for D is isomorphic 
to the orthogonal group SO(4M), and is defined by j29| ) 
with C — 1m <8> cr x ® 1. Equations ( pp| ) remain formally 
unchanged except for the replacement 2M — > AM . The 
projection U i— > C/r[/ T r _1 with r = 1m ® 1 ® icr^ takes 
the Haar measure of SO(4M) into the invariant measure 
of SO(4M)/U(2M). Ensemble averages are given by 

(f(S)) Dm = (f(UrU T r- 1 )), 



The ensemble average of a product of four U's is 

(V A 1 B 1 Uh 1 D 1 V A 2 B 2 Uc 2 D 2 ) 

4M + 1 



4M(4M- l)(4M + 2 



4M(4M- l)(4M + 2 



$Ai d $A 2 C 2 $Bi D x Sb 2 D 2 

+$a 1 c 2 Sa 2 c 1 (>b 1 d 2 Sb 2 d 1 
+Ca 1 a 2 Cc 1 c 2 Cb 1 b 2 Cd 1 d, 

Sa 1 c 1 Sa 2 c 2 Sb 1 d 2 Sb 2 d x 

+&a 1 c 2 8a 2 c 1 Sb 1 d 1 Sb 2 d 2 
+(8a 1 c 1 $a 2 c 2 + 8a 1 c 2 $a 2 c 1 )Cb 1 b 2 Cd 1 d 2 

+Ca 1 a 2 Cc 1 c 2 {$b 1 d 1 Sb 2 d 2 + Sb 1 d 2 Sb 2 d 1 ) 



The remaining calculations are the same as before. We 
obtain Sg = +1/2 for class Dili, and var(g) = 1/2 for 
class D. 



C. Conjecture for CI and Dili 

For N-systems the breaking of time-reversal symme- 
try is known to reduce var(g) by a factor of two, while 
the breaking of spin-rotation invariance causes a reduc- 
tion by a factor of four. As was said earlier, this pat- 
tern is explained by the observation that var(g) simply 
counts the number of massless modes in each universality 
class. From our experience with diagrammatic perturba- 
tion theory of the model (p j23 ) we expect the same princi- 
ple to be operative here, i.e. we expect var(<?) to be still 
determined by the number of massless modes. Indeed, 
the conductance fluctuations for the classes C and D are 
seen to be bigger than the corresponding fluctuations for 
N-systems by a factor of eight. To understand this, we 
note that there is a trivial enhancement by a factor of 
2 2 = 4 due to the transfer of two elementary charges in 
an Andreev reflection. The other factor of two can be in- 
terpreted as telling us that the number of massless modes 
a priori is twice as large: for every A-type mode, which 
is already present in the N-system, there exists an extra 
D-type (or BdG particle-hole) mode in the NS-system, 
see Sec. VI. By extrapolation we are led to the following 
conjecture: 



4 


(CI), 


2 


(C), 


1 


(Dill), 


1/2 


(D), 



(W)) 



f(U)dU. 



SO(iM) 



var(g) 



which differs from the result of Brouwer and Beenakker 
pOf who found the size of the conductance fluctuations to 
depend only weakly on whether time-reversal symmetry 
is broken or not. Note however that their result applies 
to a different situation than the one considered here (In 
their case the superconducting order parameter is homo- 
geneous in space for class CI). 
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IX. CONCLUSIONS 

In this paper we have initiated the study of a special 
family of NS-systems where the spatial variation of the 
superconducting order parameter is such that the An- 
dreev phase shift averages to zero along a typical semi- 
classical single-electron trajectory. We find such systems 
particularly interesting because the proximity effect is 
inoperative and quasiparticle states exist right at the 
chemical potential. Disorder or dynamically generated 
chaos mixes the states and leads to a novel and univer- 
sal type of level statistics within an energy window whose 
size is determined by the frequency of Andreev reflection. 
By classifying systems according to their symmetries we 
identified four universality classes, denoted by C, CI, D, 
and Dili. Time reversal is a good (broken) symmetry for 
CI and Dili (C and D), while spin is conserved (not con- 
served) for C and CI (D and Dili). For each universality 
class the joint probability distribution of the quasiparti- 
cle energy levels was given in closed form. The n-level 
correlation functions for the classes C and D were calcu- 
lated by the mapping onto a free Fermi gas on a half-line 
with Dirichlet and Neumann boundary conditions at the 
origin. The joint probability distributions of the levels for 
CI and Dili were transformed into those of the Laguerre 
Orthogonal and Laguerre Symplectic Ensembles, whose 
level statistics has been worked out completely (albeit 
with a minor computational error) by Nagao and Slevin. 

To calculate the transport properties of open systems 
in the zero-dimensional limit, we formulated a random- 
matrix model and treated it using a variant of the impu- 
rity diagram technique. An important feature we pointed 
out was the doubling of the number of low-energy modes 
in comparison with conventional normal-conducting sys- 
tems. For every A-type mode, i.e. for every BdG particle- 
particle (or hole-hole) spin-singlet or spin-triplet diffu- 
son or cooperon, there exists precisely one correspond- 
ing BdG particle-hole or D-type mode. The weak lo- 
calization correction to the average conductance for an 
NSS-geometry was calculated as a function of the "stick- 
ing probability" T and two perturbations breaking time- 
reversal symmetry and spin-rotation invariance. The 
technically more involved task of calculating the vari- 
ance of the fluctuating conductance was carried out only 
for T = 1 and the universality classes C and D, by us- 
ing an S'-matrix formalism a la Mello. We found var(g) 
to be enhanced by a factor of two relative to the rule 
var(gNs) = 4var(gisr)- We attribute this enhancement to 
the doubling of low-energy modes by the coupling to the 
superconductor. Let us emphasize that the effects we 
have studied are universal (in the ergodic limit) and are 
independent of such microscopic detail as the NS-barricr 
transmittency. 

Clearly, the present paper constitutes only a first step 
into a new and exciting research area of mesoscopic 
physics, and much more is yet to be done. Some of the 
open problems are the following, (i) We have shown how 



to solve the level statistics problem for each universal- 
ity class but more generally one might also be interested 
in the crossover between classes. Here the crossovers 
CI — > C and Dili — ► D look amenable to analytical 
techniques, since the level statistics for C and D (just 
as for the Gaussian Unitary Ensemble) maps on a free 
Fermi gas problem, (ii) Our results for the level statis- 
tics are restricted to an energy range proportional to the 
inverse mean time spent between successive Andreev re- 
flections. To access the short-time or high-energy regime 
beyond the crossover scale, our maximum-entropy ensem- 
bles need to be modified by allowing for different vari- 
ances of the random pairing and normal matrix elements, 
(iii) Although we have outlined the semiclassical inter- 
pretation of the D-type (or particle-hole) modes, a more 
detailed discussion of their role in semiclassical periodic- 
orbit theory would certainly be desirable, (iv) We need 
to extend our results for the universal conductance fluc- 
tuations to the classes CI and Dili and to arbitrary T. 
(v) While the zero-dimensional (or ergodic) limit is ad- 
equately described by the maximum-entropy ansatz, the 
diffusive and ballistic regimes necessitate a more detailed 
modelling. In particular, the nonrandom nature of the 
magnitude of the pairing field will make itself felt at short 
times. It is an open technical problem how to deal analyt- 
ically with the phase randomness of Hamiltonian matrix 
elements when their magnitude is to be kept fixed, (vi) 
We have concentrated on an NS-geometry that is partic- 
ularly easy to treat but future work will have to include 
other geometries, (vii) Last but not least, we need to 
address the nontrivial question: how large is the effect 
of residual Coulomb interactions on the D-type modes? 
There is no doubt that the short-time physics can be ad- 
equately described in a simple indcpcndcnt-quasiparticle 
picture, but at large times the coherence between parti- 
cles and holes will get cut off by Coulomb blockade and 
other correlation effects. The question is what is the time 
scale where this happens. 

We shall end on a mathematical note. According to 
Cartan, there exist eleven large families of symmetric 
spaces. Those of type II are the compact unitary, or- 
thogonal and symplectic Lie groups (A, B,C, D). The 
large families of type-I symmetric spaces are denoted by 
AI, All, Mil, DDI, CI, CII, and DHL The standard 
Wigner-Dyson universality classes derive from A (GUE), 
AI (GOE), and All (GSE), while the universality classes 
of a massless Dirac particle derive from AIII (chGUE), 
DDI (chGOE), and CII (chGSE). As we have shown, the 
new universality classes found in mesoscopic NS-systems 
exhaust the remaining large families of symmetric spaces 
except for B (the orthogonal group in odd dimensions), 
which does not occur. Thus, if we group B together with 
D, there is a bijection between the known universality 
classes of disordered single-particle systems and the large 
families of symmetric spaces. We consider this to be a 
strong indication that no other universality classes will 
be found, since any additional universality class would 
exceed Cartan's scheme and therefore would have to be 
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of a different nature. 
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APPENDIX A: CONDUCTANCE AND WARD 
IDENTITIES 

In this appendix we elaborate on the calculation of the 
weak localization correction to the conductance, Eq. (28). 
Owing to the presence of the BdG particle-hole degree of 
freedom, this calculation turns out to be much more in- 
volved than in pure N-systems. For this reason we prefer 
to control the diagrammatic expansion by means of exact 
algebraic relationships. The basic concepts used in this 
appendix have been introduced in the seminal paper [p2[ . 

To begin with, we recapitulate two Ward identities 
that will play a crucial role in what follows. Let us 
write the average retarded Green's function as (G a , a ') = 



Sa,a' (~ S^, 



4*,a'G>> where a = (fJ,,s,q) 



a composite index accounting for the site (//), spin (s) 
and particle-hole (q = p, h) degrees of freedom. (After 
averaging, the Green's function depends only on the site 
index.) The first Ward identity is immediate from the 
definitions and reads 



ft n 



-AG M (-AS M + 2*r /1 )' 



(Al) 



where G* is the average advanced Green's function, 
AG f+ = Gft-G* and AX^ = S^ — S*. A second and 
less trivial Ward-identity |52| relates the self energy S to 
the so-called irreducible two-particle vertex U: 



AS,, 



E 



U a .a' AG n 



(A2) 



The irreducible vertex U is defined as the set of all trun- 
cated four-point functions that cannot be cut by just cut- 
ting two average Green's functions, i.e. 



U 



+ 



+ 



In the following we focus on the analysis of the auxiliary 
quantity 

a' ;q'= p 

^From this the mean conductance is obtained as (cf. 



(g) = 8 ]T T^ a . 

a;q= h 



(A3) 



We start from the ansatz 



$ Q = AG, 



c + c s (-y + Cq (-y + c sq (-y + * + 



+r Al (d + 4(-) s + d g R g + 4 9 R s+g ) ■ (A4) 



(Although other expressions involving arbitrary func- 
tions of F^ seem possible, this formula does represent 
the most general starting point. Equation ( pi| ) yields 
X^T™ = 7™~ 1 X^r Al , so it is sufficient to start from 
an expression that is linear in the coefficients T^.) The 
quantity <I> satisfies Dyson's equation |32| 



(AS M - 2iY fl ) $ Q = AG P (r^ qp + J2 U a , a >$ a ^J (*), 



(A5) 



where use of the identity ( Al) has been made. To fix the 
coefficients c, ...,d, ... we subject Eq. (A5) to various 
summation procedures. For example, by taking the sum 
Ea(*) an d then using the second Ward identity (|A^), we 
obtain 



• + jd — -l/4i 



(A6) 



Seven more equations for the remaining coefficients are 
generated by performing the summations ( — 

E a (-) 9 M, EJ-) S+9 M, Ea^C*), E«r„(-)'(*), 

E Q r p (-)'(*) and Y. a r^(-) s+9 (*)- The outcome of all 
this may be cast in the form of a matrix equation 



A B 
C D 



= -2ai 



(A7) 



where c 1 = (c s ,c q ,c sq ), d? = (d,d s ,d q ,d sq ), u T = 
(0,1,0), v? = (0,0,7,0), and o x = E„r p AG M . Fortu- 
nately, it is easy to invert the 7x7 matrix appearing in 
this equation. By construction, the coefficients appear- 
ing in the subblock A(B, G, D) involve summations over 
the index a that do not (do) contain matrix elements T^. 
Since E« • • • ~ °( N )' whereas £ Q r ^ • • ■ ~ O(M), the 
coefficients appearing in A exceed those in the remain- 
ing subblocks by a large factor of O(NfM). We thus 
conclude 



+ ... 



A B 
C D 





D- 1 



0, d~D~ x v. (A8) 



The matrix D can easily be inverted as it is already of 
diagonal form. Combining Eqs. flA3|), (A6) and (A8), we 
obtain 



(g) = 16ai 



1 



27V 



(A9) 



4i 47C1 - 8i7 2 ai - U y 
where Ci = E p r M AG M AS M and 

U= ^r M AG M (-) 9 C/ Q , Q ,r^AG^(-) 9 '. 



We next decompose the self energy according to S M = 
+ <5E M into a leading order contribution plus a cor- 
rection term <5S M of 0(1/M). Anticipating that the terms 
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ST^ and U are of the same order, we obtain the prelimi- 
nary result 



A - 7 

(g) = 2MT + 2MT— '—ImST 



A(A + 7 ) 



2(A + 7 ) : 



-0(1/M), 



where ST = (M7) 1 ^2 T^ST^. We next analyze the 

building bloc ks ST and U by diagrammatic methods. 
Because Eq. (A.1C) is based on Ward identities, it au- 
tomatically incorporates the condition of unitarity. As 
a consequence, the following calculation is much simpler 
than a direct diagrammatic analysis of Eq. ( p5|) . 

To leading order in M _1 , the self-energy correction ST 
is given by 



[5] R. Gade, Nucl. Phys. B398, 499 (1993); T. Nagao and 

K. Slevin, Phys. Rev. Lett. 70, 635 (1993). 
[6] R. Oppermann, Physica A167, 301 (1990). 
[7] K. Slevin, J.L. Pichard, and P.A. Mello, submitted to 
Journal de Physique. 
(A10) I 8 ! A - F - Andreev, Zh. Eksp. Teor. Fiz. 46, 1823 (1964) [Sov. 
Phys. JETP 19, 1228 (1964)]. 



-U- 



ST 



(AH) 



where the light- (dark-)shaded region represents a D-type 
cooperon mode (a non-singular II^ -ladder). Evaluation 
of the diagrams yields (cf. the explanation in connection 
with Fig.g) 



ImST 



A(A- 7 ) 
2(A + 7 ) 



1 



MT MT + rj a 



(A12) 



The dominant contribution to the vertex correction U 
results from an A-type cooperon: 



a,ot' 



« AG M ,ry(-)« = 



4A 7 tm 



1 



MT + rjt MT + Vs 



(A13) 



By combining Eqs. flAlOf ), (|A12|) and (|A13| ), and us- 
ing the expression for the transmission coefficient T = 
4\-f /(\ + 7) 2 , we arrive at the final result given in 
Sec. IVTl. 
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